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General relativistic quasiequilibrium states of black hole-neutron star binaries are computed in 
the moving-puncture framework. We propose three conditions for determining the quasiequilibrium 
states and compare the numerical results with those obtained in the excision framework. We find that 
the results obtained in the moving-puncture framework agree with those in the excision framework 
and with those in the third post-Newtonian approximation for the cases that (i) the mass ratio of 
the binary is close to unity irrespective of the orbital separation, and (ii) the orbital separation is 
large enough (moQ < 0.02, where mo and Q, are the total mass and the orbital angular velocity, 
respectively) irrespective of the mass ratio. For mof2 > 0.03, both of the results in the moving- 
puncture and excision frameworks deviate, more or less, from those in the third post-Newtonian 
approximation. Thus the numerical results do not provide a quasicircular state, rather they seem 
to have a non-negligible eccentricity of order 0.01-0.1. We show by numerical simulation that a 
method in the moving-puncture framework can provide approximately quasicircular states in which 
the eccentricity is by a factor of ~ 2 smaller than those in quasiequilibrium given by other approaches. 

PACS numbers: 04.25.D-, 04.30.-w, 04.40.Dg 



I. INTRODUCTION 

Ground-based laser interferometric gravitational-wave 
detectors such as LIGO [J, VIRGO @, GEO600 0, 
and TAMA300 [4( are now in operation, and advanced 
detectors such as advanced LIGO will be in operation 
in the next decade and are expected to detect gravita- 
tional waves. Among many other sources, coalescing bi- 
nary compact objects such as neutron star-neutron star 
and black hole-neutron star (hereafter NS-NS and BH- 
NS, respectively) binaries are the most promising sources 
[H, Is 0] ■ This has been motivating theoretical studies of 
gravitational waves from the inspiral and merger of the 
binary compact objects. 

In the past two decades, a number of short-hard 7-ray 
bursts has been observed by the 7-ray and x-ray satel- 
lites [1, Q. However, the progenitors of these bursts are 
still highly uncertain. One of the plausible models of the 
central engine is the merger of NS-NS and/or BH-NS bi- 
naries [H [H| ■ This scenario is based on the idea that 
a system consisting of a rotating BH and hot, massive 
accretion-disk is formed as a consequence of the merger, 
and subsequently, they emit huge amount of 7-rays and 
x-rays in a short time scale. To theoretically explore 
this possibility (more specifically, to clarify the forma- 
tion process of the BH-disk system), general relativistic 
study for the merger of NS-NS and BH-NS binaries is 
probably the unique approach. This issue has been also 
motivating numerical studies for the merger of NS-NS 
and BH-NS binaries. 

In the past decade, substantial effort has been paid 
in the community of numerical relativity for clarifying 



the inspiral and merger processes of binary compact 
objects. In particular, a wide variety of simulations 
have been performed for the merger of NS-NS binaries 
[H [H, Q [H, EI ES El El, US, Eland black hole- 
black hole (hereafter BH-BH) binaries^j2p,[2l,[25|,[26, 

23, ii, JM MMMMjMj3p, HI S [Ha So, El, 

44 Ii! . |44 l45l Us 14711481149115011 . In the past two years, 



a method for computing quasie quili brium states of BH- 
NS binaries has been developed [5ll . H3, HH , and 
also, numerical simulations for the merger of BH-NS bi- 
naries have been done [H IH, [H [SO, H3 ■ However, 
these researches are still in early stage; e.g., most of the 
simulations have been performed only for a short time 
scale, and for such simulations, the inspiral and subse- 
quent merger phases are not likely to be computed very 
accurately (but see [62|, [63|). Also, it is not clear whether 
the computed quasiequilibrium states are really quasi- 
circular; i.e., it is not clear whether the eccentricity is 
sufficiently small. Numerical simulations are performed 
employing the quasiequilibrium states as the initial con- 
dition. If the quasiequilibrium state is not in a quasicircu- 
lar orbit, the results of the numerical simulation are not 
realistic. The purpose of this paper is (i) to present an 
improved study for the quasiequilibrium states and (ii) to 
compare several quasiequilibrium states obtained so far. 
In an accompany paper [63j |. we also present the latest 
accurate numerical results for the inspiral and merger of 
BH-NS binaries. 

To compute a quasiequilibrium state of BH-NS bina- 
ries, we have to employ an appropriate method by which 
the singular behavior of the BH is avoided. In most of the 
previous studies [5l|, HH, [H, HH, [55[ , the quasiequilibrium 
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state is computed in the so-called excision framework, in 
which a spherical region inside an apparent horizon is 
excised and basic equations are solved imposing plausi- 
ble boundary conditions at the excised two-sphere. In 
this paper, we employ the so-called moving-puncture ap- 
proach (see Sec. Ill All . w hich is an alternative of the ex- 
cision approach [56l . |57| . In this approach, we do not 
have to excise any region around the BH horizon nor im- 
pose boundary conditions around the BH, as shown in 
Ref. [sH • Indeed, this has been proven to be quite useful 
for computing a quasiequilibrium state [65[ and also for 
simulating binary BHs (e.g., [H 

Another possible merit in the moving- puncture ap- 
proach is that there is a flexibility for computing a 
quasiequilibrium state. In the excision method, one im- 
poses the boundary conditions at the excised two-sphere 
and at spatial infinity. The boundary conditions at the 
excised surface are usually determined by requiring that 
the two-sphere should be the Killing horizon at least ap- 
proximately. As a result, a quasiequilibrium state is com- 
pletely determined with no ambiguity, although it is not 
clear whether the obtained quasiequilibrium is really a 
quasicircular state [66|. By contrast, any boundary con- 
dition does not have to be imposed around the BH in 
the moving-puncture approach. Because of this, we have 
a remaining degree of freedom for defining a quasiequi- 
librium state. Specifically, we do not have any natural, 
physical condition for determining the center of mass of 
the system in this method. However, this degree of free- 
dom can be used to obtain a favorable, quasiequilibrium 
state. As we illustrate in Sec. lIII| quasiequilibrium states 
obtained in the excision method are not quasicircular in 
general: Namely, the eccentricity is not zero (see also 
Ref. (66|). In the moving-puncture approach, the remain- 
ing degree of freedom can be used to reduce the eccentric- 
ity, and it may be possible to obtain a quasiequilibrium 
state in which the eccentricity is smaller than that ob- 
tained by the excision method. 

This paper is organized as follows: In Sec. UH we first 
review the basic equations for computing quasiequilib- 
rium states in the moving-puncture approach. Then, 
we describe three methods for determining the center 
of mass in the moving-puncture approach. The numeri- 
cal methods for solving a quasiequilibrium state are de- 
scribed in Sec. IIIII We present our numerical results and 
compare those with other results in Sec. lIVI In Sec. lIV Cl 
we present results of numerical simulation of the inspiral 
and merger of BH-NS binaries for a particular model: We 
adopt two initial data obtained in the moving-puncture 
approach, and compare the numerical results. We then 
demonstrate that one of the moving-puncture approach is 
superior for the initial condition because the eccentricity 
is smaller than those obtain in other methods. Section 
fVl is devoted to a summary. Throughout this paper, we 
adopt the geometrical units in which G = c = 1, where 
G and c are the gravitational constant and the speed of 
light. Latin and Greek indices denote spatial and space- 
time components, respectively. 



II. FORMULATION 

In this section, we first review the basic equations for 
computing a BH-NS binary in quasiequilibrium, and then 
describe the quantities used in the analysis and the meth- 
ods for defining the quasiequilibrium state in the moving- 
puncture approach. 

A. Field equations in the moving-puncture 
framework 

To derive a quasiequilibrium state of BH-NS binaries as 
a solution of the initial value problem of general relativ- 
ity, we employ a mixture of the conformal thin-sandwich 
decomposition [qT| and conformal transverse-traceless de- 
composition of the Einstein equations. Following the pre- 
vious works, we assume that the trace part of the extrin- 
sic curvature (K =tr(-KV,-)) is zero and the three-metric 
('Jij) is conformally flat [681 ]: 7y = ip^fij, where tp is the 
conformal factor and is the flat metric. 

We define a tracefree, weighted extrinsic curvature as 

A« = ip w K ij . (1) 

Because we assume that the three-metric is conformally 
flat, this quantity is written by 

i ij = ^ + - lf j v k p*) , (2) 

where a is the lapse function, (3 k the shift vector, and V, 
the covariant derivative with respect to fy. Note that 
adding a rotational shift vector, 

0l t = (O x R) 1 , (3) 

does not change A y in the conformal flatness formal- 
ism. Here, f2 is the angular velocity vector of the binary 
and R is the coordinate vector from the center of mass 
of the binary. Thus, computation may be performed in 
any rotational frame using the same equations, by simply 
changing the boundary conditions. 

In the present formalism, the basic equations are de- 
rived from the Hamiltonian constraint, momentum con- 
straint, and the maximal slicing condition (dtK = 0) (69j 
as 

AV; = -2^ PH - \ir 7 A l0 A^ , (4) 

o 

A/3* + \^j(3 3 = Wna^f + 2A ij V j (aip- 6 ), (5) 

A$ = 2tt^ 4 ( Ph + 2S k k ) + l&Jt-tAijAV, (6) 

8 

where A = / y VjVj, $ = aip, and 

Ph = T^npnv, (7) 
f = -T^n^l (8) 
^•=T^ 7l>7j> (9) 
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Here, is the timelike hypersurface normal, Ay = 
fikfjiA kl , and the stress-energy tensor. 

For computing a quasiequilibrium state of a system 
containing BHs, we have to appropriately treat singu- 
lar behaviors of the BHs because divergent quantities 
cannot be handled in numerical computation. Most of 
the previous works for com p utin g qu asiequilibrium states 
of BH-NS binaries [H M, M, II HI have been done 
with an "excision approach," i.e., excising the region in- 
side two-sphere of apparent horizon from the computa- 
tional domain with appropriate boundary conditions at 
the excision surface. On the other hand, a "puncture" 
method [HJ was proposed by Brandt and Briigmann 
to describe multiple BHs with arbitrary linear momenta 
and spin an gula r momenta, and a "moving-puncture ap- 
proach" [24l.l28j was revealed to be quite useful in dynam- 
ical simulations. In this paper, we employ the moving- 
puncture approach, which is developed by Shibata and 
Uryu [H, [53] for the case of BH-NS binaries. In the 
puncture or moving-puncture framework we decompose 
the metric quantities into a singular part, which is writ- 
ten analytically and denotes contribution from a BH, and 
a regular part, which is obtained by numerically solving 
the basic equations. Assuming that the puncture is lo- 



cated at rp 



' p- 



we set if) and <I> as 



4> = 1 + 



$ = 1- 



M P 
2r B H 

^BH 



V, 



(10) 
(11) 



where Mp and M$ are positive constants of mass dimen- 



sion, and tbh 



is a coordinate distance from 



the puncture. Mp is an arbitrarily chosen parameter 
called the puncture mass, whereas M$ is determined by 
the condition that Arnowitt-Deser-Misner (ADM) mass 
(M ), and Komar mass agree (this condition should hold 
when the spacetime is stationary and asymptotically flat 
0, El), i.e., 



di$dS z = 



d l ^pdS l = 2-kMq. 



(12) 



Also, we decompose Ay into singular and regular parts 
as 



Aij = ViWj + VjWi - -f l3 V k W k + Kf 3 , (13) 

where Kfj is the singular part, which denotes a weighted 
extrinsic curvature associated with the linear momentum 
of the BH written by 



^. = 4- 

» 2r 2 



->BH 



k nBH 



(14) 



BH 

and l k = Xgjj/rBH- Wi denotes an auxiliary three- 
dimensional function and W l = f^Wj. Because the total 
linear momentum of the system should vanish, the lin- 
ear momentum of the BH, P„ BH , is related to that of the 



companion NS as 



P, 



BH 



(15) 



where the right-hand side denotes the (minus) linear mo- 
mentum of the NS. 

Field equations that we have to solve are summarized 
as follows: 



Acb = -2^ 5 p H --4,- 7 A lJ A^, 



(16) 



A/3 1 + -VVjP = 167r$V 3 j l + 2A iJ V i ($V~ 7 )(17) 

A V = 2^ A { PH + 2S k k ) + l^A^A^, (18) 

o 



AW, + -ViVjW j = 8W 6 Ji- 



(19) 



We note that Aij is obtained by Eq. (fl"3"|) . not by Eq. |2]), 
because A %: > is not straightforwardly defined for a = 
when we adopt Eq. In this approach, the elliptic 

equation for (3 l has to be solved because we need (3 l in 
solving hydrostatic equations (see Sec. Ill B) l. 

All the basic equations are elliptic type, and hence, 
we have to impose appropriate boundary conditions at 
spatial infinity. Because of the asymptotic flatness, the 
boundary conditions at spatial infinity r — > 00 are writ- 
ten as 



(20) 



where we assume that the equations are solved in the 
inertial frame. We note that outer boundaries are lo- 
cated at spatial infinity in our numerical computation 
(cf. Sec. IIIII) . Thus, the above condition is exactly im- 
posed. 

In contrast to the case that the excision approach is 
adopted, we do not have to impose the inner boundary 
conditions in the moving-puncture approach. This could 
be a drawback in this approach, because we cannot im- 
pose physical boundary conditions (e.g., Killing horizon 
boundary condition) for the BH. However, this could be 
also a merit, because we have a flexibility for adjusting 
a quasiequilibrium state to a desired state by using this 
degree of freedom. In Sec. Ill Dt we discuss this point in 
more detail. 



B. Hydrostatic equations 

Assuming that the NS is composed of an ideal fluid, 
we write the stress-energy tensor as 



(p + pe + p)u^u v + pg^, 



(21) 



where p is the baryon rest-mass density, e the specific 
internal energy, p the pressure, and the fluid four- 
velocity. We employ a polytropic equation of state 



p = up , 



(22) 
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where k is the polytropic constant and T the adiabatic 
index. In this paper, we set r = 2 following previous 
works [HI, [H, [53, 0, [Ell . Using the first law of thermo- 
dynamics, e is obtained as p/[(T — l)/o], and the specific 
enthalpy, h = 1 + e + p/p, is given by 



1 



„r-i 



r 



(23) 



Thus, all the thermodynamic quantities are written in 
terms of p. 

In the polytropic equation of state, all the dimensional 
quantities enter the problem only through the polytropic 
constant, and thus, are rescaled into a dimensionless form 
by normalizing in terms of the polytropic length scale, 

i? po ly^« 1/(2r - 2) . (24) 

In this paper, we present all the quantities in the dimen- 
sionless form following Refs. [HI M, El HI El . 

BH-NS binaries are never in a true equilibrium due 
to the emission of gravitational waves. However, when 
the orbital separation d is large enough, emission time 
scale of gravitational waves iow is much longer than the 
orbital period i or b as 



taw 



1.1 



( d 



V 6m 



5/2 



m 
Ap 



(25) 



where mo and p denote the total mass and reduced mass 
of the binary. Thus, except for the final inspiral phase, 
say d < 10mo, the effect of gravitational radiation re- 
action may be safely neglected, and the binary can be 
regarded to be approximately in an equilibrium state. 
Because the binaries in a close orbit should have a cir- 
cular orbit [72, l73( , the fluid configuration should be in 
hydrostatic equilibrium in the corotating frame of the 
binary system. In addition, it is believed that the mat- 
ter in most of the NS has the approximately irrotational 
velocity field for the realistic binary configurations be- 
cause the viscous time scale for the angular momentum 
transport inside the NS is much longer than the gravi- 
tational radiation reaction time scale [74], [75| . (We note 
that the actual NSs are known to have nonzero spin an- 
gular momenta and not exactly in the irrotational states. 
However, we can still approximate astrophysical NSs well 
with the irrotational velocity fields because their typical 
rotational period is lOOms-ls and are much longer than 
their typical orbital period just before mergers, ^2 — 3 
ms, and also much longer than the dynamical time scale 
of the NSs, < 1 ms.) 

The equations of relativistic hydrostatic equilibrium 
with the irrotational velocity field is derived indepen- 
dent lyin Refs. [zl, [73, u§ |. which are summarized in 
Ref. [79J. In this formulation, one assumes the presence 
of a helical Killing vector 



(26) 

For the irrotational velocity field, the relativistic vor- 
ticity is zero as 



ional velocity field, tl 
<*V = ^^{hu v ) - V„(/ra M ) = 0. 



(27) 



Using Eq. (|27p and the helical symmetric relation for the 
specific momentum £^(hu >1 ) = 0, we obtain the first in- 
tegral of the relativistic Euler equation as 



h£nU M = — C{= const). 



(28) 



To rewrite this equation into a more specific form, wc 
decompose the four-velocity in the form 



it" = u t {^ + U' 1 ), 



(29) 



where V^ 1 is a three- velocity (i.e., n^P = 0) and de- 
notes the velocity field in the comoving frame of the bi- 



nary system. Then, £ M is written as £ p 



Substituting this equation into Eq. (|28|) . we obtain [7 
h 



+ iitV 1 = C, 



(30) 



where M, = hrf i denotes the three specific momentum. 

The condition of irrotation, Eq. (|27p . implies that Ui 
is written by the gradient of a velocity potential field W 
such that 



u t = A*, 



(31) 



where Di is the covariant derivative with respect to 7^. 
Then, V 1 is written by the velocity potential as 



v l = -C 



a 1 + t^*> 



and also 



1 



h~ 2 D k ^D k m 



1/2 



(32) 



(33) 



Thus, the first integral of the Euler equation is written 
by h, 'J, and geometrical quantities. 

The equation for the velocity potential is derived from 
the continuity equation, which can be written in the pres- 
ence of the helical Killing vector as fI6\ 



0. 



(34) 



Substituting Eq. ([32f into Eq. (|34|) , an elliptic- type equa- 
tion for \E' is derived to give 



Dilpaih^D^ - u\C + I3 1 )}} = 0. 



(35) 



This equation is also written by h, VP, and geometrical 
quantities. Thus, from the first integral of the Euler 
equation and Eq. (|35[) . h and \& are computed, and sub- 
sequently, p, e, and p are determined from the equation 
of state. 



C. Global quantities 

A quasiequilibrium state is characterized by the mass 
and spin of the BH, the mass and radius of the NS, and 
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the orbital angular velocity Q. A quasiequilibrium se- 
quence should be a sequence as a function of Q with con- 
stant values of the BH mass, the BH spin, and the baryon 
rest mass of the NS, and with a fixed equation of state 
for the NS. In this paper, we assume that the BH spin is 
zero, and the BH mass is defined by the irreducible mass 
as 



M BH = 



.4 



EH 



16tt ' 



(36) 



where Aeh is the proper area of the event horizon. In 
practice, we approximate this area with that of the ap- 
parent horizon, which is computed from an integral on 
the apparent-horizon surface 



.4 



AH 



(37) 



AH 



where we use that the three- metric is conformally flat. 
In the moving-puncture approach, in contrast to the ex- 
cision approach in which the two-sphere of the apparent 
horizon is readily known to be the excision surface, the 
apparent horizon has to be determined by a numerical 
computation. For finding the apparent horizon, we use 
the algorithm developed by Lin and Novak (see Ref. [8(| 
for details). 

The baryon rest mass of the NS is defined by 



-gdV, 



(38) 



where g is the determinant of the spacetime metric g^ v . 
In this paper, we always present the mass in polytropic 
units as 



M NS 

m s - Mb 



M, 



poly 



(39) 



whereby Mg S is the baryon rest mass for the polytropic 
constant k = 1. 

In the following, we often compare numerical results 
with those derived by the third post-Newtonian (3PN) 
approximation for two points masses [8ll . |82| . In such a 
case, we have to define the total mass mo and the reduced 
mass n of the system (i.e., we have to define each mass 
of the binary component). For the BH mass, we use the 
irreducible mass. For the NS mass, we use the ADM 
mass of an isolated NS Mj^ M with the same baryon rest 
mass. Thus, 



M 



BH 



M NS 



7M-BH M NS 



m 



(40) 
(41) 



Note that for a nonspinning BH, M£ r is equal to the 
ADM mass of the BH in isolation. 

The ADM mass of the whole system Mq is defined by 



Mo = ~ (f d^dS 1 



(42) 



and the Komar mass is 



Mi 



1 

4tt 



d i adS i 



(43) 



-J- / - d^)dS\ (44) 



Komar — 



Note that equating these two masses results in Eq. (fT 
We also define the binding energy of the binary as 

E b = M - mo. (45) 

The ADM linear momentum of the system is 

1 



Pi 



8tt 



K l3 dS\ 



(46) 



where we assume the maximal slicing K = 0. This is set 
to be zero in the present work. The angular momentum 
of the system around the center of mass of the binary 
may be defined by 

Ji = ^-Ctffc / (X j K kl - X k K jl )dSi, (47) 

where X % is the coordinate vector from the center of 
mass. 



D. Free parameters 

To compute a quasiequilibrium configuration for a 
given equation of state, we have to fix free parameters. 
Each configuration is determined when we fix (i) the 
baryon rest mass of the NS, Mg , (ii) the mass ratio, 
Q, or equivalently, the irreducible mass of the BH, M™, 
and (iii) the separation between the BH and the NS, 
d. Then, the other parameters such as Mp, M$, and 
P^ 11 are automatically determined in the computation: 
The puncture mass Mp is arranged to give a desired irre- 
ducible mass, the mass parameter Mj, is determined by 
the condition that the ADM mass and the Komar mass 
agree (see Eq. (|12]>). and the linear momentum of the BH, 
P™, is determined by the condition that the total linear 
momentum of the system should vanish []see Eq. (fl5|)]. 

There are also free parameters associated with the con- 
figuration of the NS; the integration constant C, which 
appears in the first integral of the Euler Eq. (|30"|) , and the 
angular velocity of the binary Q. These arc determined 
by fixing the configuration of the binary and the baryon 
rest mass of the NS. Specifically, we fix the rest mass of 
the NS to determine C and fix the location of the center 
of the NS to determine VI. Here, the center of the NS is 
defined as the position at which the following condition 
is satisfied [Ellzl: 



a In ft 



dX 



= 0, 



(48) 



(Xns.Ins.O) 



where Ans and Ins are the distances in the x and y 
directions from the center of the NS to the rotational 
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axis, respectively. In the actual calculation, we arrange 
the value of the specific enthalpy at the center of the NS, 
h c , instead of C, since it is easier to implement. The 
value of C is determined by h c and the values of metric 
quantities at the center of the NS. 

The final remaining task in the moving-puncture 
framework is to determine the center of mass of the sys- 
tem. The issue in this framework is that we do not have 
any natural, physical condition for determining it. (By 
contrast, the condition is automatically derived in the 
excision framework [5l|, HH, [EH, [H, [55j] , although it is not 
clear whether the resulting quasicquilibrium is a quasi- 
circular state [66|.) In our previous paper [56|, [57]], we 
employed a condition that the dipole part of ip at spa- 
tial infinity is zero (hereafter we refer to this condition 
as "dipole condition"). However, we found that in this 
condition, the angular momentum derived for a close or- 
bit of fimo > 0.03 is by ~ 2% smaller than that derived 
by the 3PN approximation [8l|, [13] for Q = 3. Because 
the 3PN approximation should be an excellent approxi- 
mation of general relativity for a fairly distant orbit as 
Slrn-o « 0.03, the obtained initial data deviates from the 
true quasicircular state, and hence, the initial orbit would 
be eccentric. 

In the subsequent work [58|, we adopted a condition 
that the azimuthal component of the shift vector /3 V at 
the location of the puncture (r = rp) is equal to — fi; 
i.e., we imposed a corotating gauge condition at the lo- 
cation of the puncture. In the following, we refer to this 
condition as the U (3 V condition." In this paper, we first 
present the numerical results in the latter condition be- 
cause it is a physical condition and gives a slightly better 
result than the dipole condition does. 

As shown in Sec. HVi however, the angular momentum 
derived for a close orbit of fimo > 0.03 in this method 
is still by ~ 2% smaller than that derived by the 3PN 
relation for a large mass ratio Q > 2. The disagreement 
is larger for the larger mass ratio. Such initial conditions 
are likely to deviate from the true quasicircular state and 
hence the orbital eccentricity is large (e.g. Sec. IV C of 
Ref. [6l| for numerical evolution of such initial data). 
This also suggests that the /3 V condition is not suitable 
for deriving realistic quasicircular states. 

In this paper, we further propose a new condition in 
which the center of mass is determined in a phenomeno- 
logical manner: We impose the condition that the total 
angular momentum of the system for a given value of 
f2mo agrees with that derived by the 3PN approxima- 
tion. This can be achieved by appropriately choosing the 
position of the center of mass. With this method, the 
drawback in the previous two methods (i.e., the angular 
momentum becomes smaller than the expected value) is 
overcome. We refer to this method as "3PN-J condition" 
in the following. 



III. NUMERICAL METHODS 

In this section, we summarize our numerical scheme 
for computing a quasiequilibrium state. 

A. Methods of computation 

Our numerical code is based on the spectral meth- 
ods library LORENE [831 ] . In the spectral methods, any 
quantity is denoted by the expansion into a complete set 
of polynomials. The feature in the spectral methods is 
that the error of this expansion decreases exponentially 
as the number of the used complete set of polynomials 
increases, at least for continuous functions. Furthermore, 
irregular functions like the rest-mass density of the NS, 
for which the spatial derivative jumps at its surface and 
is not straightforward to expand into a complete set of 
polynomials (known as "Gibbs phenomenon"), can be 
treated appropriately by employing a multidomain spec- 
tral method. We use two sets of spherical-like computa- 
tional domains, each of which is centered on each object. 
The one for which the domain center is located at the 
puncture is composed of one nucleus, which is sphere at 
the center, several shells, and the external domain, which 
extends to spatial infinity. The other has almost the same 
structure, except that the outer boundary of the nucleus 
is deformed to fit the surface of the NS. With this com- 
putational domain, irregular profiles of the density are 
contained only in the domain boundary, and hence, no 
Gibbs phenomenon arises, if the density decreases suf- 
ficiently smooth at the boundary, e.g., for the T = 2 
polytropic equation of state. We also locate the outer 
boundaries at spatial infinity by employing a radial coor- 
dinate, which is obtained by a transformation u = 1/r in 
the outermost domain. With this treatment, the exact 
boundary conditions at spatial infinity can be imposed. 
Because the multidomain method is employed and field 
equations are solved for many domains, we split the field 
equations into the "BH part" and the "NS part" like the 
way described in Appendix. A of Ref. [H| . 

Our computational domains centered on the BH are 
divided into 8 domains and each of them is covered by 
N r xNgX — 41 x 33 x 32 collocation points. Similarly, 
the domains for the NS are divided into 6 domains and 
each of them is covered by N r x Ng x = 25 x 17 x 16 
collocation points. Here, N r , N$ and are the num- 
ber of collocation points for radial, polar, and azimuthal 
directions, respectively. 

B. Iteration procedure 

Numerical solutions of quasiequilibrium states are ob- 
tained by iteratively solving the basic equations described 
in Sec. [TTJ Here, we briefly summarize our procedure of 
the iteration. As described in Sec. Ill Dl the calculation 
should be performed to give correct Mg S , M™, and d. 
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First of all, we need to prepare an initial trial solution 
for the iterative procedure. For this purpose, we superim- 
pose a Schwarzschild solution in the isotropic coordinates 
and a spherical NS. Then, our iterative procedure is as 
follows: 

1. Determine the orbital angular velocity f2 by using 
Eq. dUD- 

2. Determine the location of the center of mass (rota- 
tional axis) of the system. Because the coordinate 
separation between two objects are initially given, 
the positions of both objects relative to the center 
of mass are also determined. 



3. Solve the equations described in Sees . Ill Al and HTbI 
in each domain. 

4. Adjust Mp to fix the irreducible mass of the BH to 
a desirable value. After this procedure, determine 
M$ so that the condition (fT2l) is satisfied. 



5. Adjust the maximum enthalpy of the NS, h c , at the 
center of the NS, to fix the baryon rest mass of the 
NS. 

We repeat this procedure until a sufficient convergence is 
achieved. As a measure of the convergence, we monitor 
the relative difference between the enthalpy field of suc- 
cessive steps. Typically, we stop the iteration when the 
following condition is satisfied: 



J2\ 1 ~ h t^/h7jl\<io- 



(49) 



where n denotes the iteration step, and (i,j,k) the col- 
location points of (r, 9, (p). 



IV. NUMERICAL RESULTS 

Throughout this paper, we characterize a quasiequi- 
librium sequence by two parameters; the baryon rest 
mass of the NS, Mq , and the ratio of the irreducible 
mass of the BH to the ADM mass of the NS in isola- 
tion, Q = M^ r H /M^f] M . We focus on the sequences 
of M£ s = 0.14, 0.15, and 0.16 following Ref. [54[. For 
these cases, the compactness of the NSs is C — 0.1321, 
0.1452, and 0.1600, respectively, implying that we choose 
moderately large values for the compactness. Note that 
the maximum value of M^ s is about 0.18 for the T = 2 
polytropic equation of state. At that value of Afg S , the 
compactness is about C = 0.21. Here, the compactness 
is defined by 



C 



"ADM,0 
Rq 



(50) 



where Rq is the circumferential radius of the NS in iso- 
lation. For the mass ratio, we choose 1 < Q < 5. 



A. Binding energy and total angular momentum in 
the f3 v condition 

Figure [1] plots the binding energy and total angular 
momentum as functions of Qtuq for Q = 1, 3, 5, and 
Mg S = 0.15 in the (3 V condition. For comparison, the 
results in the excision approach [54| and in the 3PN ap- 
proximation are also plotted. 

For large orbital separations (small values of Q,m^ < 
0.02), the results obtained in the moving-puncture 
method (with the (3 V condition) agree well with those 
derived by the 3PN approximation and in the excision 
approach [5~ij irrespective of the mass ratio. For smaller 
separations, however, the degree of agreement among 
three results depends on the mass ratio. For Q = 1, 
three results agree within ~ 1% error. By contrast, for 
Q = 3 and 5, a significant deviation of order ~ 10% 
arises among three results for Vim® > 0.03. In par- 
ticular, for Q = 5, the results in the moving- puncture 
approach (with the (3 V condition) disagree significantly 
with other two results (see lower panels of Fig. [J). Be- 
cause the tidal effect does not play an important role 
and the orbital velocity is at most ~ 0.1c for a fairly dis- 
tant orbit of Qmo ~ 0.03, the numerical results should 
agree with the results in the 3PN approximation, at least 
approximately. This implies that the quasiequilibrium 
state obtained in the moving-puncture approach (with 
the f3 v condition) is not in a quasicircular orbit for the 
close orbit of £lm > 0.03 and for Q > 3; it would be 
an eccentric orbit. Also, quasiequilibrium states in this 
moving-puncture approach appear to be inferior to those 
in the excision approach in that they show systematic 
deviations from other two results. 

One possible reason for this deviation may stem from 
the condition for determining the center of mass of the 
system. This point is explored in detail in Sec. lIVBl An- 
other possible reason is that the BH might have nonzero 
spin in the present approach. Indeed, in the excision 
framework, it has been found that the "leading-order ap- 
proximation" leads to a slightly spinning BH [84, 85]. To 
obtain a strictly nonspinning BH, a computation has to 
be performed with an improved method for determining 
the spin angular velocity of the BH. Motivated by this 
fact, we measured the spin of the BH using a method 
proposed by Cook and Whiting [86j. However, we find 
that the BH has negligible spin of order S/M? T < 10~ 5 , 
and therefore, the deviation between the results in the 
moving-puncture method (with the /3 V condition) and 
others is not caused by the spin of the BH. 

Numerical results of binding energy and total angular 
momentum in the moving-puncture method as a function 
of Otoo for different compactness of the NS are plotted 
in Fig. O This shows that the feature of the results de- 
scribed above holds irrespective of the compactness of 
the NS. For smaller values of the compactness (i.e., for 
smaller values of Mg S ), the binding energy and the an- 
gular momentum are slightly larger for a given value of 
Qttiq > 0.04. This is due to the fact that for less compact 
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FIG. 1: Left panels: Binding energy Eb/mo as a function of fimo for Mj S = 0.15 and Q = 1 (upper panel), 3 (middle panel), 
and 5 (lower panel). The solid (red) and dashed (green) curves show the results obtained in the moving-puncture method 
with the I3 V condition and in the excision method [5J], respectively. The dotted (blue) curve denotes the result in the 3PN 
approximation [8l| . Right panels: The same as the left panels but for the total angular momentum J/mo as a function of flmo. 



NSs, the tidal-deformation effect on the NS plays an im- 
portant role in increasing these quantities in close orbits 
[87| . We note that, in the moving-puncture approach, it 
is possible to obtain the sequences of quasiequilibria for 
the NSs with the compactness C < 0.18 for the T = 2 
polytropic equation of state. Meanwhile, our computa- 
tions do not show adequate convergence for the binaries 
containing more compact NSs, because such NSs are close 
to the most compact ones allowed by the given equation 
of state, C ~ 0.21 in this case, and are difficult to com- 



pute accurately. For the T > 2 equation of state for 
which the maximum compactness is larger than 0.21, we 
are able to obtain NSs of compactness ~ 0.2. 

Before closing this subsection, we point out the fol- 
lowing issue found from Fig. [TJ The angular momentum 
for a given value of Slmo in the numerical results is al- 
ways smaller than that in the 3PN approximation for 
Qirio > 0.03. This holds for the results not only in the 
moving-puncture approach (with the f3 v condition) but 
also in the excision approach. For a fairly distant orbit of 
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FIG. 2: The same as Fig. [T]but for M^ s = 0.14, 0.15, and 0.16, and for Q = 3 in the moving-puncture method. 



flmo ~ 0.03, the 3PN approximation should be an excel- 
lent approximation of general relativity, and hence, pro- 
vides a highly accurate result. This implies that angular 
momentum in the numerical results is smaller than that 
for the real quasicircular state. As shown in Sec. IIV CI 
this is indeed the case. Figure[T]indicates that for numeri- 
cal simulation, these quasiequilibria may not be good ini- 
tial conditions, and an improved quasiequilibrium would 
be favorable as the initial condition of numerical simula- 
tion. 



B. Effect of the center of mass 

As mentioned in Sec. IIID| we have no definitive guid- 
ance for determining the position of the center of mass in 
the moving-puncture framework. In the excision frame- 
work, the position is automatically determined so that 
the ADM linear momenta should vanish. By contrast, 
in the moving-puncture framework, this condition was 
already used to determine another free parameter, Pj BH . 

As described in Sec. Ill D| we have at least three meth- 
ods for determining the center of mass, and numerical 
results depend strongly on them. In this subsection, we 
compare these numerical results. 

In Fig.[3l we show sequences of Mg S = 0.15 for Q = 1, 
3, and 5 obtained by three different methods for the cen- 
ter of mass, "shift," "dipole," and "set to 3PN" denote 
the results derived in the f3 v , dipole, and 3PN-J condi- 
tions, respectively. 

The numerical results in the (3 V and dipole conditions 
show a similar behavior. For Q = 1, these results agree 
approximately with those in the 3PN approximation, and 
for Q = 3 and 5, the deviation from the 3PN results be- 
comes significant as pointed out in Sec. IIV Al The devia- 
tion from the 3PN results is slightly smaller for the results 
in the [3 V condition for Q — 3 and 5 than those in the 
dipole condition. This indicates that the (3 V condition is 
slightly better for computing unequal-mass BH-NS bina- 
ries in quasiequilibrium, although the deviation from the 



3PN results is larger than 1% for fimo > 0.03. 

The sequence obtained in the 3PN-J condition shows 
rather different behavior. By definition, the angular mo- 
mentum agrees with the results obtained in the 3PN ap- 
proximation. As a result, however, the binding energy 
becomes larger than that in the 3PN approximation, in 
contrast to the results in other two methods. A possible 
interpretation for the excess of the binding energy is that 
junk gravitational radiation or nonstationary kinetic en- 
ergy (e.g., oscillation of a nonstationary BH) are included 
in the data. However, if these nonstationary components 
are radiated away during numerical evolution, this initial 
data could provide an approximate quasicircular state. 



C. Assessing quality of quasiequilibrium by 
numerical simulation 

The results presented in the previous subsection show 
that the quasiequilibrium states computed in three con- 
ditions of the moving-puncture approach are not in qua- 
sicircular orbits for Q ^ 1, but rather are likely to be 
in eccentric orbits. However, if the eccentricity is small 
enough, it quickly approaches zero during evolution, re- 
sulting in that a realistic binary, i.e., a (n approximately) 
quasicircular state, is provided. To assess circularity of 
quasiequilibrium states as the initial condition of nu- 
merical relativity, we performed a numerical simulation, 
choosing a binary of M^ s = 0.15 and Q = 3 and 5. 
In this subsection, we present the results of the numeri- 
cal simulation. In this experiment, we adopt the initial 
data obtained in the j3 v and 3PN-J conditions of the 
moving-puncture approach. For all the data, the angular 
velocity is chosen to be fimo rs 0.033. The numerical sim- 
ulation was performed using our code SACRA, in which 
an adaptive mesh refinement algorithm is implemented. 
The formulation, gauge condition, and numerical scheme 
adopted in SACRA are the same as described in Ref. [6l[ . 

Figure 2] plots evolution of the orbital angular velocity 
as a function of retarded time. Here, the angular velocity 
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FIG. 3: Left panels: Binding energy E/mo as a function of fimo for M^ s — 0.15 and Q = 1, 3, and 5 in the moving-puncture 
approach with three different conditions for determining the center of mass of system. The dotted curve denotes the result in 
the 3PN approximation. Right panels: The same as the left panels but for the total angular momentum J/rriQ as a function of 



is calculated by the quadrupole mode of gravitational 
waveforms by 



1 \^ i (l = m = 2)\ 



dt^ 4 (l = m= 2) 



(51) 



where ^4 denotes the outgoing part of the Newman- 
Penrose quantity. This figure illustrates that the orbit 
of the binaries is eccentric for all the cases, reflecting 
the fact that the circular orbit is not provided initially. 



For the initial data obtained in the (3 V condition, the 
eccentricity appears to be of order 0.1. Moreover, the 
eccentricity does not reduce to zero even at the onset of 
merger. We note that for these cases, the binary spends 
in the inspiral phase for 4-5 orbits before the onset of the 
merger. Nevertheless, the eccentricity is not sufficiently 
reduced by the gravitation radiation reaction and a qua- 
sicircular orbit is not achieved before the merger. The 
reason for this is that the initial eccentricity is too large. 

By contrast, for the initial data obtained in the 3PN-J 
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for C — 0.15. trot = approximately denotes the merger time. The solid and dashed curves denote the results in the 3PN-J 
condition and the f3 v condition, respectively. The dot-dotted curve is the result in the 3PN approximation (Taylor T4 formula; 
e.g., Ref. SI). 



condition, the eccentricity appears to be much smaller 
than those for other two cases. This suggests that the 
initial eccentricity is smaller. 

By comparing the numerical results with that obtained 
in the 3PN approximation (dot-dotted curve), we find 
that the modulation amplitude of the angular velocity 
is A(Qmrj) ~ 0.004-0.006 even just before the merger 
(at fimo ~ 0.05) in the /3 V condition. Eccentricity of 
the orbit is approximately estimated to be 2AS1/317 for 
a slightly eccentric orbit. Thus, the eccentricity for these 
cases is - 0.05-0.08. 

By contrast, the modulation amplitude is at most 
A(Omo) ~ 0.003 for the initial data given in the 3PN-J 
condition. Furthermore, the modulation is suppressed to 
be < 0.001 just before the onset of merger (Qmo > 0.05) 
in this case, indicating that the eccentricity is < 0.01. 
This shows that the moving-puncture approach with the 
3PN-J condition is superior for providing the initial data 
for the numerical-relativity simulation. 



D. Mass-shedding limit 

As analyzed in detail in Refs. [H, |54|, the NS is 
strongly subject to tidal deformation by a companion 
BH outside their innermost stable circular orbit (ISCO), 
if the mass ratio of the system is small enough or the ra- 
dius of the NS is large enough. In the case that the tidal 
field of the BH is strong enough, the Lagrangian point 
enters inside the surface of the NS. Then the NS cannot 
be in equilibrium any longer because mass shedding oc- 
curs from the inner edge of the NS's surface. We here 
determine the condition for the onset of mass shedding 
in the moving-puncture approach with the (3 V condition. 
As shown in the previous two subsections, the quasiequi- 
librium states obtained in this approach are not quasi- 
circular states but slightly eccentric ones. This seems 
to be also the case for the quasiequilibrium states ob- 
tained in the excision approach employed in Refs. [53l.l54j|. 
Even though there exists nonzero eccentricity in the data, 
we think that it still deserves to reinvestigate the mass- 



shedding limit in the moving-puncture approach and to 
compare the results with those obtained in the excision 
approach. 

In the spectral methods one cannot determine any 
equilibrium configuration for a star with irregular sur- 
face shape (i.e., with a cusp). At the onset of mass shed- 
ding, the inner edge of the NS has a cusp, and hence, 
it is not possible to accurately compute quasiequilibrium 
states of the NS in close orbits with a BH. Thus, we infer 
the orbital separation (or angular velocity) at the onset 
of mass shedding from quasiequilibrium states of slightly 
more distant orbits than the mass-shedding configuration 
has. For this procedure, we define a mass-shedding indi- 
cator x [Z§| as the fraction of the radial derivative of the 
enthalpy at the NS surface, 



X 



_ ( (d(\nh)/dr] 



cq 



\(d(lah)/dr) po i e 



(52) 



where subscripts "eq" and "pole" imply that the partial 
derivative is taken with respect to the radial direction 
connecting the centers of the NS and BH in the equa- 
torial plane and along the polar direction, respectively. 
For the infinite separation, the NS becomes spherical and 
X = 1 , whereas x decreases with increasing the degree of 
tidal deformation and eventually becomes zero at the on- 
set of mass shedding. Because our code does not converge 
to give quasiequilibrium for x ~~ * 0, we first compute a 
sequence of x for close orbits, and then, derive a fitting 
formula for the curve of x as a function of orbital sepa- 
ration. By using this fitting formula, we determine the 
orbital separation for x = 0. 

Figure[5]plots the mass-shedding indicator x as a func- 
tion of fimo for the sequences of Mg S = 0.15, and Q = 3 
and 5. For comparison, we also plot the results obtained 
in the excision method. This figure shows that the curves 
obtained in the moving-puncture approach agree approx- 
imately with those in the excision approach. This indi- 
cates that the mass-shedding limit would be determined 
with a small error even in the presence of spurious eccen- 
tricity of order ~ 0.1. 

In the works performed in the excision approach 
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FIG. 5: The mass-shedding parameter \ versus Qmo for the 
sequences of M^ s = 0.15, and Q = 3 and 5. We also show 
the results obtained in the excision method [541 ] . 

[Hll |54| . the ISCO is also determined by finding the ex- 
tremum of the binding energy (or the total angular mo- 
mentum). They find that for M^ s = 0.15 and Q = 5, 
the binary reaches the ISCO before the mass-shedding 
limit is reached, whereas for Q — 3, the mass shedding 
occurs before the binary reaches the ISCO. In the present 
moving-puncture approach, we do not find any extremum 
even for Q — 5. This is primarily due to the fact that 
this approach is not suitable for computing quasicircular 
states for very close orbits. 

We also compare our results of mass-shedding limit 
with that of the dynamical simulations [63j. Since it is 
difficult to determine the onset of mass shedding in the 
dynamical simulations, we determine the value of Qmo at 
the time when the matter in the NSs is first swallowed by 
the BHs. For the case of Q = 3, the matter is swallowed 
by the BH at ttmo ~ 0.09, which is larger than the mass- 
shedding limit obtained in this paper and in other works 
[53], [HI]. However, two results are consistent because the 
mass shedding occurs before the matter falls into the BH; 
flrriQ ~ 0.09 shows the upper limit for the value of l]mo 
at the mass shedding. On the other hand, for the case 
of Q = 5, the matter is first swallowed by the BH at 
flrriQ ~ 0.1. This simply implies that the binary reaches 
the ISCO at ftm - 0.1. 



V. SUMMARY 

We numerically derive new general relativistic 
quasiequilibrium states of BH-NS binaries in the moving- 
puncture approach. Basic equations for gravita- 
tional fields are solved in the mixture of conformal 
thin-sandwich decomposition and conformal transverse- 
traceless decomposition, and hydrostatic equations are 
solved under the assumption of irrotational velocity field 
and employing the polytropic equation of state. 

In the moving-puncture approach, no definitive physi- 



cal condition is present for determining the center of mass 
of the system. We propose three conditions for defining 
the center of mass and compare the resulting quasiequi- 
librium states with the 3PN results and results obtained 
in the excision approach. 

Sequences of quasiequilibrium states are computed for 
the mass ratio 1 < Q < 5, and for M^ 3 = 0.14, 0.15, 
and 0.16 (each of which corresponds to the compactness 
C = 0.1321, 0.1452, and 0.1600). For large orbital separa- 
tion with Omo < 0.02, the results in the moving-puncture 
approach agree with the results obtained in the 3PN ap- 
proximation and in the excision approach, irrespective 
of the conditions for determining the center of mass and 
irrespective of the mass ratio of the binary. Thus, such 
quasiequilibrium states are likely to be (at least approx- 
imately) quasicircular orbits for which the eccentricity 
is approximately zero. By contrast, for small orbital 
separation with Qtuq > 0.03, the results in the moving- 
puncture approach with (3 V and dipole conditions for de- 
termining the center of mass systematically deviate from 
the 3PN results, in particular, for large mass ratio. The 
angular momentum for the resulting quasiequilibrium is 
always smaller than that in the 3PN results, and hence, 
the quasiequilibrium appears to contain a nonzero eccen- 
tricity. As an alternative method for determining the 
center of mass, we propose a method in which it is de- 
termined by requiring that the angular momentum for a 
give value of the angular velocity should agree with that 
in the 3PN approximation. 

To assess the circularity of the quasiequilibria com- 
puted in different approaches, we performed numerical 
simulation for a chosen model; Mg S = 0.15, Q = 3 and 
5, and Qmo w 0.033. The numerical simulation was per- 
formed for two quasiequilibrium states prepared in the 
(3 V and 3PN-J conditions. We find that the quasiequi- 
libria computed in the j3 v condition have eccentricity of 
~ 0.05-0.08. During the simulation, the eccentricity re- 
duces due to gravitational radiation reaction, but in ~ 4- 
5 orbits, it does not reduce to < 0.01 even at the onset 
of merger. By contrast, the eccentricity of the quasiequi- 
librium state computed in the 3PN-J condition is by a 
factor of ~ 2 smaller. For such case, the eccentricity re- 
duces approximately to ~ 0.01 in ~ 4-5 orbits, and the 
resulting eccentricity at the onset of merger appears to be 
< 0.01. Therefore, with the quasiequilibrium prepared in 
the 3PN-J condition, it is feasible to compute a realistic 
gravitational waveform at least from the final a few inspi- 
ral orbits to the merger phase (see an accompany paper 
[63[ for detailed numerical results). 

In this paper, we study BH-NS binaries in which the 
BH is not spinning. If the BH has a substantial spin, the 
quasiequilibrium state of BH-NS binaries will be modified 
by the spin-orbit interaction effect (e.g. Ref. |82j). We 
plan to study the effect of the BH spin in the next work. 
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